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Abstract 

We study the numerical resolution of the time-dependent Gross-Pitaevskii equation, a non-linear 
Schrodinger equation used to simulate the dynamics of Bose-Einstein condensates. Considering 
condensates trapped in harmonic potentials, we present an efficient algorithm by making use of a 
spectral Galerkin method, using a basis set of harmonic oscillator functions, and the Gauss-Hermite 
quadrature. We apply this algorithm to the simulation of condensate breathing and scissors modes. 
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I. INTRODUCTION 

The experimental realization of Bose-Einstein condensation has prompted much 

work on the study of the dynamics of these condensates. From the theoretical side, many 
interesting results have been obtained using the Gross-Pitaevskii equation (GPE) ^0,0], 
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with the normalization condition ||\l/(t)||^2 = 1 Vt, to describe the order parameter \1/ (also 
called the condensate wave function) of condensed bosons of mass m, interacting via a 
contact potential described by the scattering length a, and eventually confined by an external 
potential V^xt- Even though the Gross-Pitaevskii equation is based on the approximation 
that all bosons are in the condensed phase (T = K), direct comparison between theoretical 
and experimental results have shown that, in many cases, solutions of the GPE contain the 
essential physics of the underlying phenomena 

m 

\m. This non-linear Schrodinger 
equation (NLSE) has been used, in its time-dependent form, to investi gate many aspects of 
the dynamics of Bose-condensed gas, such as the formation of vortices IllH . the interference 



between condensates 
a few. 



12j . of the possibility of creating atom lasers 



to mention only 



Most of theses and other numerical studies of the time- dependent GPE are based on grid 
methods, i.e., discretize the spatial coordinates on a grid of points, the resulting differential 
equation b eing usually solved by Crank-Nicholson or split-operator Fourier methods (see. 



e.g., Refs. 
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llS^). We must point out that, while much care must be taken in 
solving Eq. (Q) because of the non-linearity, we find, to our dismay, that many authors give 
results calculated with the time- dependent GPE without even specifying what method they 
have used for their numerical simulation. 

In this article, we wish to focus our attention on the case where the Bose-Einstein con- 
densate is in a (possibly anisotropic) harmonic trap, i.e.. 



VUX. Y, Z) = -m {ulX' + + ^'^') (2) 

The method we propose is based 
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which is the case for most experimental set-ups 
on the spectral decomposition of \1/ on a basis of harmonic-oscillator wave functions. In such 
a representation, the kinetic -|- trapping potential part of the Hamiltonian is diagonal. The 
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non-linear part is computed by forward and backward transformations from the spectral 
to a grid representation. By judicious use of the Gauss-Hermite quadrature, this can lead 
to an algorithm that is more efficient than those based on grid methods. Although this is 
akin to DVR methods based on Hermite polynomials, which have been successfully used for 
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our method is distinct, since our 



the time-independent and time-dependent GPE 
Hamiltonian is expressed in the spectral representation for both the kinetic and potential 
operators. 

We expose in Sec. HT] our spectral method and the resulting algorithm. We then present 
different time evolution schemes that can be used in combination with the spectral method. 
We finally give in Sec. II VI some results that can be obtained from the numerical simulation 
of the time-dependent GPE, namely the study of condensate breathing and scissors modes. 

II. SPACE DISCRETIZATION 

To simplify the calculation, we will first rescale Eq. ^ in the three spatial dimensions 
{X, Y, Z) and in time, 

1/2 

] 



X = { ) X (3a) 



mujy 



Y = [-—] y (3b) 



t = —T. (3d) 



We also introduce a new wave function ip defined as 

^{t,X,Y,Z)=A^{r,x,y,z) 
and, considering the normalization condition 

[ \^{t,X,Y,Z)fdXdYdZ = l\ft, 



we choose 
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such that 



(r, X, y, z) 1^ dxdydz = 1. 



The Gross-Pitaevskii equation therefore becomes 



dip 



^ (4) 



with 

1/2 



X = AnaNr^^] . (5) 

Coordinate z should be chosen such that uJz is the greatest of the three frequencies [this is 
related to the arbitrary choice of the scaling factor in Eq. pd]) ]. 

As all the physical parameters have been absorbed in the non-linear parameter A, cal- 
culations with the same A can correspond to results for different species, but in diverse 
experimental conditions. We can define acceptable lower and upper bounds for A by consid- 
ering the effective range of the different physical parameters. Considering only cases where 
the interparticle interaction is repulsive, i.e., a > and therefore A > 0, at the lower end we 
can consider a small "^He* condensate (m = 4.0 a.m.u., a = 302 a.u. [23]) of = 10^ atoms 
in a highly anisotropic uo^uOy/ujz = 27r x 10^^ Hz trap, giving A ~ 1.3, while for a bigger 
N ~ 10^ condensate of heavy atoms such as '^'^Rb (m = 86.9 a.m.u., a = 106 a.u. [24']), A 
can reach 10^ for isotropic traps. In the following, we will restrict our study to A in the range 
1-10'^, considering that the Thomas- Fermi approximation can be used for greater values of 

A H- 

A. The spectral Galerkin method in ID 

For pedagogical purposes, we first explain our numerical method on the simple case of 
the one-dimensional NLSE 

t^{t,x) = Ho^{t,x) + X\ij{t,x)\^ij{t,x) (6) 



with 

^0 = + ^x^- 



Id' 1 2 



2 9x2 2 

Extensions to the three-dimensional case will be detailed in the next section. 
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Denoting by ijj{t) the function x t-^ ilj{t,x), it can be proven j25| that if 



[ ^ < +00, f x'^\x{x)\'^dx < +oo\- 
Jr (jx Jr J 



Eq. (jni) with initial condition ipo has a unique solution in C°([0, +oo[, X)nC^{[0, +oo[, 
and that both the norm 

1/2 



L2 



and the energy 



\ilj{t,x)\ dx 



X 



E = {Hom,m) + '-^ I mt,x)\'dx 



are conserved by the dynamics. A variational formulation of Eq. (jH}, supplemented by the 
initial condition ip(t = 0) = ipo where i/jq G X, reads 



Search E C%[0,T], X) n C\[0,T], L\R)) such that 



(7) 



Numerical solutions can then be obtained by approximating problem ((7j) with a Galerkin 
method: a finite dimensional subspace Xn of the infinite dimensional vector space X being 
given, we consider 

Search ■ipN' ^ C^{[0,T], Xj^) such that 

dt 

Denoting by (00, • • • , 4'n) an orthonormal basis of Xjv for the scalar product and by 
C{t) = (c„(t))o<ri<Ar the vector of C^^-*^ collecting the coefficients of ipNif) in the basis 
(00, • • • ,07v), i.e., 

N 

i'N{t,x) = y^C„(t)0n( 



n=0 



problem (jH]) can be reformulated as 



Search C G C\[0,T],C^^^) such that 
dC 

i^{t) = hC{t) + \F{C{t)) 
C(0) = Co 



(9) 



where Co are the coefficients of ipo and h the matrix of Hq in the basis (0o, • • • , (Pn) 



[Co]n = (^0, 0n)L2, hnm = (-f^O^m, 4>n) , 

and where the function F is defined by 

N 

nc)n= E (10) 

with 

■^klmn 

The efficiency of a direct implementation [2a, |22| of the Galerkin method described above 
is very poor: the calculation of the integrals Ikimn (which can be precomputed if the basis 
is small enough that the integrals can be stored in memory) scales as 0{N'^Np) where Np 
is the number of grid points of the quadrature method, and the computation cost for one 
evaluation of the function F scales as A^^ [for each of the coefficients, 0{N^) operations 
are needed]. 

Our aim is to show that the Galerkin method becomes very efficient if (00, . . . , ^at) are 
the + 1 lowest eigenmodes of the harmonic oscillator Hq. In this case, indeed, the vector 
F{C) can be computed exactly (up to round-off errors) in 0{N'^) operations. Let us recall 
that the eigenmodes (0n)neN of Hq read 

where Ti.n{x) is the n-th Hermite polynomial [28], and that they satisfy 

ifo0n = En(j)n, with En = n+^. 

In such a basis, the matrix h is therefore diagonal: h = Diag(£'o, . . . , En)- In addition, for 
any C G C^^"^, one has 

F{C)n= [ mx)\^^{x)Mx)dx, (11) 

where ip{x) = "^^=0 (^n<Pn{x). The key point is now that for any n < N the integrand in (jllll 
is of the form Q{x) e~^^^ , where Q{x) is a polynomial of degree lower or equal to AN; each 
of the + 1 integrals can therefore be computed exactly with a Gauss-Hermite quadrature 
formula involving 2A^ Gauss points 2^- More precisely, we have, for any polynomial Q of 
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degree lower or equal to AN, 

/ Q{x)e~''''dx = ^ WkQ{xk) 
k=i 

where {xhl are the roots of the Hermite polynomial I-L2N+1 and where {wk} are convenient 
weights j30|. By a change of variable in integral (lllj). it follows that 

Spectral Galerkin methods are usually not very efficient jsil; but they can be in the specific 
case of the NLSE we are interested in because of the special form of the nonlinearity. 

Let us now denote by P G Ai{N + 1, 2N + 1) the matrix collecting the values of the basis 
functions {4>n)o<n<N at the Gauss points {xk)i<k<2N+i- 

p 

-' nk 



and by Wk = Wke^^ / An efficient algorithm for the computation of F{C) for a given 
C e C^+i reads: 

1. Compute the vector \E' G C^^"''^ defined by 

^ = • C. 

2. Compute the vector S G C^^^^ coefficient by coefficient along formula 

3. Compute 

F{C) =P-E. 

The vectors C and \E' are the representation of the wave function in the spectral basis 
{'Pn}o<n<N ^^^1 space (at the 2N + 1 Gauss points {x^/v^}), respectively. Steps 1 

and 3 of the above algorithm scale quadratically in (these are matrix- vector products), 
and step 2 scales linearly in A^. We therefore end up with an algorithmic complexity in 
0{N^). 

In practice, the function C ^ F{C) is called one or several times at each time step; of 
course, the matrix P as well as the weights Wk can be precomputed once and for all and 
stored in memory. 



B. The spectral-Galerkin method in 3D 



Let us now turn to the 3D setting and consider the rescaled equation 

i){t,x,y,z) + A|^/'(t,x,y, z)\^'ilj{t,x,y,z) 



(12) 



with 



Hn(x) 



1 1 2 , , 1 92 1 2 , , 
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For A > 0, a global- in-time existence and uniqueness result is available for Eq. (jl2j) with 
initial condition ilj{t = 0) = ipQ and 

i^oeX = {xeL\m'), Vxe{L\m'))\ {x' + y' + z'Y^\ e L\m')] . 

On the other hand, it is well known that finite-time blow-up may be observed for A < and 
for some initial conditions As stated above, we focus here on the case where A > 0. 

Following the same lines as in the Sec. Ill Al the approximated wave function tpN{t) is 
expended on the spectral tensor basis set 



One therefore has 



Ny 

^N{t,X,y,z)= X] C„^n,n,(i)0n,(a;)0„,(l/)0„,(2). 



(13) 



nx=0 ny=0 nz=0 



The equation satisfied by the three index tensor C = [cn^nyn^] in the Galerkin approximation 
formally has the same expression as in ID, 

dC 

t^{t) = hC{t) + \F{C{t)), 
the linear operator h now being defined by 



[hC]nxnynz En^nyn^Cn 



with 



X ThyTlz 



2 / cuz 
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and the non-linear function F{C) by 



{x, y, z)\^'iIj{x, y, z)(f)n^ (a;)</)„ {y)4>n, {z)dxdydz, 



where ip{x, y, z) is given by Eq. (fT5|) . 

Let us denote by {xk}i<k<2N^+i-> iVk} i<k<2Ny+i-> i<k<2Nz+i roots of the Hermite 



polynomials H2n^+i, 'H2Ny+i, ^2Af,+i and {wl}^ 



<k<2N^+l^ 



<fc<2Af„+l' 



Mil 



<k<2N:,+l 



the 



associated summation weights. Let us also introduce the matrices G M.{Nx + 1, 2A^^ + 1), 
Py e M{Ny + 1, 2Ny + 1), eM{N, + 1, 2N, + 1) defined by 

[Px]n^k^ = (j)nAxkjV2), [Py]nyky = (pUyiVky / ^2) , [Pz]n,k, = KX^kJ^), 

and the weights 



2 

wl e '-y 



w 



The following algorithm for the computation of F{C) scales in 0{NNxNyNz) where N 

max{N,,Ny,N,): 

1. Set ^^^^ = C 



2. Compute vl>ff«,^ 



z\nzkj^ nxnyUz 



n^=0 

Ny 

3. Compute ^l^^^ = J^lPyU^y^'J 

ny=0 

4. Compute ^f^k. = E t^J"^'^^^". 



SRR 



5. Compute E^l\ = wi wf wl l^f^^. ?^kTk 



6. Compute Ef^ff. 

7. Compute SgJ^^ 



8. Compute EfJ^^^ 



fcz=l 

2Ny + l 

E\ p 1 ■=:-Ri?S' 
l^ylnyky^k^kyTlz 

2N^+1 



0{N^NyN^) operations 
0{N^N^Nz) operations 

0{NlNyNz) operations 

0{NxNyNz) operations 
0{N^NyN^) operations 

0{N^N^Nz) operations 

0{NlNyNz) operations 



9. Set F(C) = H^^^. 

In the above formulation, the superscripts 5* and R stand for spectral and real space repre- 
sentations respectively. In other words, steps 2-4 constitute the successive transform of the 
wave function from the spectral basis to a spatial representation on the series of points of 



9 



the Gauss-Hermite quadrature. The non-hnear term of the Hamiltonian is then calculated 
in this spatial representation (step 5), while steps 6-8 correspond to the inverse transform 
back to the spectral basis. It is this procedure of forward and backward transformation that 
allows us to obtain a much better scaling than the implementation of Eq. (jlUj) . 

The scaling of the above algorithm [0{N'^) if = Ny = N^] has to be compared with 
the scaling of FFT based algorithms which scale in 0{Np\og{Np)) where Np is the number 
of grid points per direction. The main interest of the spectral method is that for a similar 
accuracy, the number of spectral basis functions per direction (here denoted by A^) can 
usually be chosen much smaller than the number Np of grid points per direction. This is 
especially true when the problem considered displays a symmetry in one or more of the 
directions, in which case the basis set used in the Galerkin approximation Eq. (jl3|) can be 
restricted to even harmonic oscillator functions (in the corresponding direction). We will 
come back on this important feature of the spectral method in Section IIVI 



C. Exploiting spherical or cylindrical symmetry 



When 



ujz the one-particle Hamiltonian possesses spherical symmetry. If the 



initial condition ipQ = ip{t = 0) has the same symmetry, then the wave function ip{t) is 
spherical symmetric for any t > 0: ■ip{t, x, y, z) = ijj{t, r) where r = (x^ + + z'^Y^'^, is the 
radial coordinate. Eq. (j3]) leads to the effective ID dynamics 



dijj 



Let us now define the function 



1 d 



d 



27r ril)(t, r) if r > 

— \/2tt rip{t, —r) if r < 0. 



It is easy to check that x actually satisfies the ID NLSE 



HoX + A 



dt """^ ' ■'27rr2 



X- 



(14) 



Besides, for any t > the function x(t) : r x(t,r) is odd and belongs to L^(]R) since 



\x{t,r)\'^dr= I 4:7rr^\ip{t,r)\'^dr = 1. 

'0 
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It can thus be expanded on the odd modes of the harmonic oscillator: 



+00 



X{t,r) = ^Cn(i)02n+l(r). 



n=0 



A spectral Galerkin approximation can now be used. The vector C(t) e C^+^ collecting the 
coefficients {ck{t))o<k<N of the approximated wave function 

N 



XN{t,r) = ^Cn{t)(j)2n+l{r) 



n=0 



obeys once again a dynamics of the form 

dC 



Here 



and 



i—{t)^hC{t) + XF{C{t)). 



h = Diag(£;2n+i), with E2n+i = 2n + -, 



[F{C)]n = / 



■X{r)(f>2n+i{r)dr, 



where x{r) = Zl„=o '^"<^2n+i(?^)- As for any < n < A^, 02„+i(r) = rP2n{r)C' where 
is a polynomial of degree equal to 2n, it follows that the above integrals can be computed 
exactly with 4A" Gauss points. 

Let us now turn to the cylindrical symmetry when (for instance) cUx — ojy and when the 
initial data reads ipo{x,y,z) — i/jo{r,z) with r — {x'^ + y^)^/^. In this case, the cylindrical 
symmetry is preserved by the dynamics so that for any t > 0, i/j(t,x,y,z) = ip{t,r,z) and 
the time evolution of ip(t,r, z) is then governed by the 2D equation 



1 d 



2r dr V dr 



d 



1 92 



2dz^ 



(15) 



set on the spatial domain 1R+ x IR. Defining a new function x(t, r, z) on the space domain 
by 

tp(t, r,z) if r > 

ip{t, —r, z) if r < 



it occurs that x satisfies 



1 92 



2 9r2 



1 d' 



2 9^2 



1 d 



ijjz 2r dr 



X 



(16) 
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on the space domain M^, and that, by construction, the function r i— r, 2;) is even. A 
spectral Galerkin approximation is obtained by expanding the wave function on the spectral 
tensor basis set 

The coefficients {cn^nJo<nr<Nr,o<nz<Nz of the expansion are solution of an equation of the 
same form as above, 

^^it) = hC{t) + XF{C{t)). 
The main difference is that in this case, the linear map h takes into account the operator 

_J_d_. 

2r dr ' 




Let us remark that the scalar product , j is well defined since the first derivative 

of 4>2nr is of the form rP2n,.{r)e^'^^^'^ where P2nr is a polynomial of degree 2nr; in addition, 
it can be computed exactly by numerical integration with 2^^ Gauss points. It is worth 
pointing out that the "Hamiltonian" in (fT?)|) is not self-adjoint because of the term — 
and that the norm of x{t) is not a conserved quantity; on the other hand, the norm 
of x{t) for the measure r dr dz is conserved. 



III. TIME DISCRETIZATION 



When a spectral Galerkin method is used to discretize the space variables, one ends up 
with a finite dimensional dynamical system of the form 

dC 

z—{t) = hC{t) + XF{C{t)), (17) 
with initial condition C(t = 0) = Cq. We then use a basic fourth-order Runge-Kutta 



method 32| to solve Eqs. ()17p. Let us mention that, as the Hamiltonian character of the 
NLSE is preserved by the spectral Galerkin discretization, it would be possible to resort 
to symplectic methods 0]; such algorithms, which are particularly advised for long time 
evolution, are however not tested in the present work. 

We will also use a grid method, based on the split-operator method, to serve as a bench- 
mark for the spectral algorithm we have just detailed. We recall below the main features of 
this approach. 
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The wave function at time r + Ar can be obtained from the wave function at r according 



to 



ip{T + Ar) =U{t,t + AT)i;{T), (18) 
with the propagator U{t, t + Ar) being expressed, for sufficiently small intervals Ar, as 



f/(r, r + Ar) = exp [-ii^ (r) Ar] 



(19) 



where H{t) is the Hamiltonian of Eq. (Q. As the potential and non-linear components of 
the Gross-Pitaevskii Hamiltonian do not commute with the kinetic operator, we apply the 
split-operator method js^ to obtain 



exp [— ■j/7(r)Ar] = exp 



~iT- 



,At 



exp [-i {V + Xl-^f) At] exp 



-iT — 
2 



+ 0(Ar=^), (20) 



with T the kinetic operator and V the trapping potential. The middle term is diagonal 
in position space, while the kinetic part is diagonal in momentum space. A Fast Fourier 
Transform is thus used before application of the kinetic operator, followed by the inverse 
transform. Note that if the intermediate wave function at time r + Ar is not needed, the 
.wo .ucce.s,ve M„etic op.ato. half-steps can be eo.bi.ed. F... a p.evous study Q. 
it appears that the split-operator method is the fastest algorithm for solving a NLSE on a 
grid. 



IV. RESULTS 



The first test we perform is the propagation of the ground stationary state (obtained 
from the time-independent GPE solved by a method based on the Optimal Damping 
Algorithm while monitoring the value of the coefficients c(r) of the ex- 

pansion (fT!?|) . For the spherically symmetric case, we require that the relative error 
on the Co coefficient (which has the largest absolute value) be inferior to 10^*^, i.e., 
)|2| / |co(r = 0)1^ < 10-^ Vr G [0,100]. This criterion also results in 



Co{r - Co(r 



an absolute error of all coefficients |c„(r) 



< 10 ^. We have also checked that 



2 2 

the phase of the coefficients is correct, by calculating |c„(r) — c„(r = 0)e~*'^'n / |c„(r)| , 
where fx is the chemical potential of the ground stationary state of the GPE j6|, and this 
value indeed is less than 10~^^. 
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In this ID case, we need = 20 basis functions for A = 100, and the resulting time-step 
for the Runge-Kutta propagator is Ar = 0.005. If A = 1000, the basis set used should be 
slightly larger, = 26, with a smaller time step Ar = 0.0025 to insure that the above 
error criteria are met. The resulting propagation time up to r = 100 is 8.9 s for A = 100 
(calculated on an Athlon 1.2 GHz processor running under Linux, using the NAG Fortran 
90 compiler at the -0 level of optimization) and 28.3 s for A = 1000. If we double the size of 
the basis set, we get a CPU time of 32.9 s for A = 100, showing the expected 0{N^) scaling 
of the algorithm in ID. 

Comparing now with the grid method described in Sec. IIIH we use Np = 64 grid points 
in the range —8 < r < 8. The time step used is Ar = 0.00025, resulting in a propagation 
time of 10.3 s, which is slightly longer than what we obtain using the Runge-Kutta method. 

We now apply our algorithm to study the dynamics of trapped condensates. Referring 
again to the spherically symmetric case, we start with the stationary ground state for an 
isotropic trap frequency u. We then let this initial state ipo evolve in a trap of frequency 
uj/2, as illustrated in Fig. corresponding to an experiment where the frequency of the 
potential trapping the condensate would be instantaneously reduced by a factor of 2. The 
corresponding time-evolving wave function \ilj{t,r)f is shown in Fig.|21 for A = 10. We must 
note that the values of A we give correspond to the condensate in the initial cu-frequency 
trap, the effective value being used for the time evolution is thus scaled by 1/V2 [see Eq. (jSJ], 
while r is rescaled with respect to the final trap frequency uj/2. We can see the "breathing" 
of the condensate as it expands and recontracts in the trap. 

It is also interesting to look at the effect of the value of the non-linear parameter on 
the breathing frequency of the condensate, as seen in Fig. El First, we note that the initial 
density at the center of the trap is lower for bigger values of A, which is expected because 
of the corresponding higher interparticle repulsion. Starting from an unperturbed harmonic 
oscillator (A = 0), for which the complete cycle time is r = 47r with recurrences every r = vr, 
we observe that the oscillation frequency of the condensate in the trap increases with a 
greater value of A. 



For the 3D case, we will study the scissors mode j39|, |4^ of a trapped condensate. We 
consider a pancake-shaped condensate, formed in an anisotropic trap with = Uy u^, 
see Fig. |3] The y and z axes of the trap are instantaneously rotated, at t = 0, by angle 
6 around the x axis. The condensate then starts to oscillate in the trap, leading to the 
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FIG. 1: Trapping potentials to (solid line) and uj/2 (dashed line) used to simulate the breathing 
modes of a condensate. The wave function |V'('")|^ of the stationary state for potential uj with 
A = 100 is also given (dotted line). 

so-called scissors mode. 

Using the parameters of the experiment of Marago et al. [4^, we first determine the 
stationary state for a condensate of = 10^ ^''Rb atoms in a trap with ujz = 255 Hz, 
^x/^z = ^y/^z = resulting in a value A = 147.1. The condensate is then tilted by an 

angle of 6* = 3.6°, with the trapping frequency reduced by 2%, resulting in a new value 
of A = 148.6. We then calculate the free evolution of this tilted condensate. We report, in 
Fig. 13 the angle between the condensate (as determined by the main inertia axis) and the y 
RXIS, clS cl function of time for the free evolution of the condensate in a trap. The oscillation 
frequency, in these conditions, is found to be 1.105 (in rescaled units), corresponding to 
276 Hz. 

This simulation was done using a basis set of = 29 functions in each dimension, using a 
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FIG. 2: "Breathing" of the condensate after expansion from a trap of frequency uj to lv/2. The 
density profile is given as a function of time r (scaled with respect to the final trap frequency 

Lij/2) for an initial A = 10 ground stationary state. 

time step Ar = 0.005. The calculation time for a propagation of duration r is then 1735 s. 
The main advantage of using a spectral Galerkin method, as noted in Sec. Ill B| is that in this 
case we can restrict the basis set in the x dimension by using only even harmonic oscillator 
functions, since the reflection symmetry with respect to the yOz plane is conserved. The 
number of functions is thus reduced to A^^^ = 15, resulting in a decrease of CPU time to 
~ 1030 s for 1 r. This compares favorably with the grid method, for which an equivalent 
calculation with 64 x 64 x 64 grid points takes ~ 1700 s (using the same time step and grid 
spacing as for the ID grid method). 
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FIG. 3: "Breathing" of the condensate after expansion from a trap of frequency a; to a;/2. The 
value of the wave function in the center of the trap, \'^{r = 0)|^, is given as a function of time r 
(scaled with respect to the final trap frequency io/2) for A equal to (solid line), 10 (dotted line), 
100 (dashed line), and 1000 (dot-dashed line). 

V. CONCLUSION 

We have presented the application of a spectral Galerkin method to the numerical so- 
lution of the Gross-Pitaevskii equation, describing a Bose-Einstein condensate trapped in 
a harmonic potential well. This method is based on the decomposition of the condensate 
wave function on the a basis set of eigenmodes of the harmonic oscillator, while the nonlinear 
term in the GPE is calculated using the Gauss-Hermite quadrature. The resulting algorithm 
scales in O(iV^) for a full 3D problem (where N is the number of basis functions used per 
direction), which is slightly worse than the O(A'^logA^p) scaling obtained for grid-based 
Fourier methods. Nevertheless, the required number of basis functions needed for a given 
problem can be much smaller than the number of grid points Np, allowing for fast and effi- 
cient calculations using the spectral method. We have shown how the propagation in time 
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FIG. 4: Representation of the study a condensate's scissors mode. The condensate is initially tilted 
with respect to the trap's y and z axes by an angle Q. 

can be carried out using a Runge-Kutta method on a set of coupled ordinary differential 
equations. 



This method is akin to the DVR approach 
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which is relies on the fact that matrix 



elements of the nonlinear term can be approximately evaluated to sufficiently high accuracy 
using an iV-point rule based Gauss quadratures. Our approach has the advantage that, for 
the basis chosen, there are no approximations in the computation of these integrals. Let 
us however remark that the same property can hold within the DVR method by a suitable 
choice of weighted polynomials. The main distinction between the usual DVR approach and 
our method is that we treat the kinetic and potential terms of the Hamiltonian conjointly, 
as detailed in Sec. HH 

We have successfully applied our algorithm to simulate two different dynamical aspects 
of trapped BECs. Making use of the spherical symmetry of an isotropic trapping potential, 
we used an effective ID equation to study the breathing of a condensate that is allowed to 
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FIG. 5: Time evolution of the angle 6 for the scissors mode. The crosses correspond to the angle 
resulting from the time-dependent calculation, along with the corresponding Hi 9 = 3.6 cos (l.lOSr) 
(dashed line). 

expand from more confining trap to a looser one. In the 3D case, we have looked at the 
scissors modes of a pancake-shaped condensate, for which the trapping potential is suddenly 
rotated along one axis. 

Future work will focus on the implementation of better time-evolution algorithms on 
our spectral method and on its possible parallelization. Extensions will also be made to 
consider other terms in the Gross-Pitaevskii Hamiltonian, such as the potential created by 
the interaction with a laser field, or coupled Gross-Pitaevskii equations used in the simulation 



of two-species condensates |41j or of the formation of molecules in atomic condensates 
43|. 
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